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We revisit, both numerically and analytically, the finite-time blowup 
of the infinite-energy solution of 3D Euler equations of stagnation- 
point-type introduced by Gibbon et al. (1999). By employing the 
method of mapping to regular systems, presented in Bustamante 
(2011) and extended to the symmetry-plane case by Mulungye et 
al. (2015), we establish a curious property of this solution that was 
not observed in early studies: before but near singularity time, the 
blowup goes from a fast transient to a slower regime that is well 
resolved spectrally, even at mid-resolutions of 512^. This late-time 
regime has an atypical spectrum: it is Gaussian rather than expo¬ 
nential in the wavenumbers. The analyticity-strip width decays to 
zero in a finite time, albeit so slowly that it remains well above the 
collocation-point scale for all simulation times t < T* — 10“®°'^'^, 
where T* is the singularity time. Reaching such a proximity to sin¬ 
gularity time is not possible in the original temporal variable, because 
floating point double precision (r; 10“^®) creates a ‘machine-epsilon’ 
barrier. Due to this limitation on the original independent variable, 
the mapped variables now provide an improved assessment of the 
relevant blowup quantities, crucially with acceptable accuracy at an 
unprecedented closeness to the singularity time: T* — t m 10“^^®. 
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pare various numerical methods, for the purposes of accurately 
solving the system and diagnosing blowup. In this regard we 
reinvestigate the stagnation-point-type solutions (with infinite 
energy) of the 3D Euler equations found to have analytic 
solutions exhibiting hnite-time blowup [7], now employing a 
novel new method that maps nonlinearly the time and fields 
to a globally regular system [S]. By doing this we uncover 
curious late-time behaviour of the Fourier spectrum such that 
in the new variables the solution remains well spectrally con¬ 
verged until far beyond the time at which the original system 
could hope to achieve due to floating point precision proximity 
to the singularity time T*. We hnd that the mapped variables 
maintain acceptable levels of error in the main blowup quan¬ 
tities such as the and L® norms of the vortex stretching 
rate at this extreme closeness to T*. We begin by formulating 
both the original and mapped equations and reviewing and 
updating the analytic and asymptotic results for blowup. We 
then present a thorough investigation of the Fourier spectra 
of the solution, followed by error analysis of our numerics and 
an assessment of singularity time and proximity to it. 


Abbreviations: 3D, three-dimensional; ODE, ordinary differential equation; GPU, 
graphics processing unit; BKM, Beale-Kato-Majda 

T he open question of regularity of the fluid dynamical 
equations is considered one of the most fundamental chal¬ 
lenges of mathematics and physics [T]. While the viscous 
Navier-Stokes equations have more physical relevance, the in¬ 
viscid Euler equations present the greatest challenge and ex¬ 
hibit the most extreme behaviours. For this reason the numer¬ 
ical study of possible hnite-time blowup is typically concerned 
with these inviscid equations. The three-dimensional Euler 
equations for an incompressible huid of unit mass density with 
velocity held u(a;, y, z, t) G R®, in a time interval t G [0, T), can 
be expressed as: 


Formulation 

We consider a class of exact solutions of the 3D Euler equa¬ 
tions presented by Gibbon et al. [S]. Writing \i{x,y,z,t) — 
{uxix,y,t),Uy{x,y,t),z'y{x,y,t)) we obtain 


^+u,.V,7 

= 2(7^)-7^ 

[2 

duj _ 

^ + Uh • VhCc; 
ot 

= 70 ;, 

[3 


where Uh{x,y,t) = {ux{x,y,t),Uy{x,y,t)) denotes the “hori¬ 
zontal” component of the velocity held at the symmetry plane 
(2 = 0 ), Vh = {dx,dy) denotes the “horizontal” gradient op¬ 
erator, uj is the vorticity dehned as 


-I- u ■ Vu = —Vp, V • u = 0. [1] 

Periodic boundary conditions are commonly assumed in a fun¬ 
damental domain [0, 27r]® . Beale et. al. provided a blowup 
criterion based on the maximum vorticity (inhnity norm) pro¬ 
viding an important reference for diagnosing singularities nu¬ 
merically (BKM theorem). The literature surrounding the 
numerical assessment of hnite-time blowup in 3D Euler is ex¬ 
tensive and will not be reviewed in detail here (see [3] , |4] and 
for references pertinent to the current work see the introduc¬ 
tion in M), sufhce to say that a fundamental difficulty of this 
important problem is the lack of analytic solutions or any a 
priori knowledge of asymptotic behaviour. A secondary obsta¬ 
cle is that the spatial collapse associated with intense vortex 
stretching results in numerical solutions becoming unresolved 
beyond a certain time (e.g., loss of spectral convergence). It 
is therefore these authors’ view that it is imperative to de¬ 
vise a framework with nontrivial blowup dynamics and where 
analytic solutions are known in order to validate and com¬ 


- edxUy OyUx^ 


7 is the stretching-rate of vorticiW which using the incom¬ 
pressibility condition in equation Q can be dehned as: 

'y{x,y,t) = -\7h-Uh{x,y,t), [4] 


and 


(/(*, P, t)) = y^ f{x,y,t)dxdy 


denotes the spatial average over the periodic 2D domain. 

Constantin [7] solved for 7 along characteristics (and for 
vorticity lj, which can be found a posteriori), proving that the 
stretching rate 7 would blow up in a hnite time, with explicit 
formulae for the singularity time which conhrmed the accuracy 
of the numerical blowup predictions in [9]. A BKM m type 
of theorem was established by Gibbon |10| where the blowup 
time T* is dehned as the smallest time at which 



||7(-,t')IUdt' = 00 , 
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where || 7 (-,i)||oo is the supremum norm of the vorticity 
stretching rate. 

Bustamante [S] and later Mulungye et. al. introduced 
the following ‘mapped’ fields and ‘mapped’ time: 


7map(a;,y,T) = 




ll7(-,t)l|oc 

uj{x,y,t) 


"W = / ll7(',i')lloodt', 
^0 


ll7(-,t)l|oc’ 

This transformation is bijective for t < T*. The mapped fields 
satisfy the following PDE system: 


dr 

dUmai 

dr 


“h U-niap ' V'ymap = 2('^niap) 'Tmap 

“hTooO'map{l 

T ^map ' ^^map — O^map ^map 

"K^oo ^map {1-2(7 Lp)} 


where CToo = sign 7 (X.y(f),f) is the sign of 7 at the position 
X-,(f) of maximum | 7 (x, t)\. The initial conditions used in this 
study are: 

7 o(a;, y) = uJoix, y) = sm(a:) sm{y) . [ 8 ] 


Analytical solution of the stagnation-point-type 3D Euler 
flows. System S-O is an exact solution of 3D Euler equa¬ 
tions (albeit with iimnite energy), as derived originally by Gib¬ 
bon et. al. [6]. Ohkitani and Gibbon [9] performed a numer¬ 
ical study at resolution 256^, supported with simulations at 
resolution 1024^, which provided evidence of a finite-time sin¬ 
gularity at f ~ 1.4. Higher resolution was not needed due to 
the fact that spectral convergence was observed during most 
of the simulation time. 

Constantin introduced a method for finding analytically 
the blowup quantities (e.g. ||7(-,t)||oo, ( 7 ^)) and established 
that there is a finite-time singularity. While it is possible to 
obtain the asymptotic behaviour of the blowup quantities us¬ 
ing the notation in [7], we will discuss this in the context of 
our method for the purposes of simplicity of presentation. 
The equation of motion for S{t), the function upon which pro¬ 
vides the solution for all fields along characteristics, (see for 
more details) is, following spatial integration over the initial 
condition: 


S = 


4 [K{S^)]^ ’ 


S{0) = 0, 


[9] 


where K{y?) is the complete elliptic function of the first kind 
such that S can be interpreted in this case as the modulus 
of the elliptic function. The singularity occurs when S = S*, 
where S* = —1/ inf 70 = 1. Thus, S(t) goes from 0 at f = 0 to 
1 at f = T*, where T* is the singularity time defined as 


A Z' ^ 

T* = — / [A'(5'^)]^dS« 1.418002734923858875062234. 
Jo 

The singularity is dominated by the infimum of 7 (x, t) and 
for these initial conditions we can identify || 7 ||oo = — inf 7 , so 
Goa = — 1 will be used throughout. 

Solving the above ODE for S{t) is in principle feasible nu¬ 
merically to any desired accuracy, although we do not know 
if a closed form is available. However, it is possible to obtain 
some exact formulae for || 7 (- , t)||cxD and r(f) in implicit form 
through S{t)-. 


- inf 7(x,t) 

xeT2 


ll7(-,i)IU = 


TT^ [{S+l)K{S^) 

AS {1 - S'^ ) K {S^f 

[ 11 ] 


( 7 ^) = 


TT^ [K{S'^) - E{S^)] [2E{S^) - (1 - S'2) K{S^)] 

3252 (1 - 


Notice that it is possible to find S{t) and 


[ 12 ] 

ll 7 (-,i(T))l|oo as 
functions of r via first inverting equation (10 1 to obtain S{t), 
and then using this to obtain || 7 (-,f(T))||tx). In numerical im¬ 
plementations at values of r greater than about 33, this re¬ 
quires the use of arbitrary-precision computations, provided 
by commercial packages such as Mathematica as S{t) becomes 
to within double floating point precision to 1 . 

Defining Z = — In(l^) , the following formulae are valid 
asymptotically as S ^ 1 : 


T* -t 


ll7(-,i)l|oo 

( 7 ^) 


^{Z^ + 2Z + 2) , 

8 y z3 J ’ 

TT^e^^ (Z-2\ 

128 V / ’ 


[13] 

[14] 

[15] 


Also, we obtain r ~ Z — In (^) to lowest order. This latter 
formula can be inverted in the asymptotic region of interest, 
giving 



where W-i is a branch of the Lambert (or Product Log) func¬ 
tion. By combining the above it is possible to obtain explicit 
asymptotic expressions for T* — t and || 7 (- , t)||oo in terms of 
the mapped time r. 

These asymptotic formulae are very useful in practice. At 
r = 5 the above asymptotic formula for f(r) has a relative er¬ 
ror of about 10“® and for || 7 (- , t)||oo, a relative error of 10 
By r > 20 the asymptotic formulae above are accurate to 
double precision (10~^®). 


Numerical solution of original and mapped systems 

We solve the evolution equations for both systems numeri¬ 
cally using a standard pseudospectral method implemented on 
GPUs using CUDA [5]. Dealiasing is carried out using Hou’s 
exponential filter exp (—36 (2fe/A)^®) (for a given spatial reso¬ 
lution N) [TT] and a fourth-order Runge-Kutta scheme solves 
in time. Adaptive time-stepping, {dt = (iT/|| 7 (-, f)||cK)), is used 
for the original equations and uniform steps of dr are used in 
the mapped system with the resulting distribution of temporal 
data roughly equivalent. 

Spectra and analyticity strip. To investigate the spatial col¬ 
lapse associated with the singularity in 7 , we consider a de¬ 
tailed analysis of the one-dimensional spectra of stretching 
rate 7 constructed from spherical shells: 

E{k,t)= l7(k,t)r 

fc— i<|k|<fe+i 


T{t) 




Our first observation of the evolution of the spectrum is that 
there are two timescales in evidence. An initial burst can 
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be observed with a flux towards intermediate k which is re¬ 
distributed across the modes. Provided N > 256 this initial 
phase remains well resolved and lasts only until r « 25. There¬ 
after there is a slow cascade from small k. In fact in original 
variables the initial phase is until T* — t ~ 10“^°. As will be 
shown, this is too early to establish certain asymptotic trends. 

It was also found that, due to the lack of direct energy cas¬ 
cade to large k, an accumulation of round-off error propagates 
up-scale. The result is a small quantity of spurious energy be¬ 
tween the large scales and the truncation wavenumber. The 
amount of this spurious energy is resolution dependent lead¬ 
ing to an ill-converged spectrum. This issue was remedied 
by applying a small amount of hyperviscosity on the large 
wavenumbers, namely adding the term 

/i^2, for|kj> 200 [16] 

to the right hand side of the Fourier transform of equation 
f6f and Q. Numerically a Crank-Nicolson scheme was used 
on this term for stability. Figure shows the profile of the 
spectra at t = 5, 10 and 25 for N = 1024 and 4096, each 
with u = 10“® and u = 0. This demonstrates that the hy¬ 
perviscosity gives a well converged spectrum while leaving the 
large scale modes unaffected. The error in the bulk quantity 
(7map) is unchanged (figure not shown), however applying hy¬ 
perviscosity to all modes leads to a significant error increase. 
Interestingly the late time profile does not have the typical 
shape we might expect | 12l or that which is assumed previ¬ 
ously in this system [5], namely 

E{k,t) = [17] 


In fact, as can be seen in figure]^ (bottom right at late times), 
the profile assumes a more Gaussian shape, 


E{k,t) = [18] 


This late time spatial form has been missed in previous work 
[9] on this system as it only arises after the initial burst, which 
does have the shape, and persists to sufficiently close to 

T* to render it next to inaccessible without the mapped vari¬ 
ables. To ensure the convergence of the initial burst phase we 
first perform a le ast-squares fitting procedure to the spectrum 
with ansatz ( |17[ ). Figure shows the fitted di for the early 
burst phase. The plot shows two resolutions (N = 1024 and 
2048) which are essentially indistinguishable. From Figure!^ 
the exponential part of the profile of E{k, t) is preceded (in ^ 
by the Gaussian shape. The downscale flux associated with 
the slackening of the exponential part (tfi decreasing) estab¬ 
lishes the Gaussian profile in its wake. The result is that, 
while the trend in (5i suggests an exponential decay in r (at 
early times), in reality the ansatz ( |17[ ) ceases to be a valid 
analyticity measure due to the addition of a large scale Gaus¬ 
sian spectrum. The cross over regime is indicated by negative 
values of for 13 < r < 23. 

To a naly se the true late time behaviour (t > 25) we fit with 
ansatz (18l. Figure [t shows the behaviour of 82 as a fun ction 
of r. Strikingly, the decay is now very slow (see Eq. (22|). 

Using a classical method, it is possible to obtain a rigor¬ 
ous upper bound for the supremum norm of stretching rate in 



T 



T 

Fig 1. Time evolution of the analyticity distance of the Fourier spectrum 
of . Top: profile at early times based on equation jl7| , showing the initial burst. 
Bottom; late time profile based on equation jl8| showing the slow Gaussian cascade 
along with the estimate (5 ^ yj' k jr, which saturates inequality |22|. 
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r 

Fig. 2. Time evolution of the exponent, n(t(T)), and constant factor, C*(t(T)) 
(normalised to coincide with the mapped variables), of the fit ansatz equation | |18| of 
the Fourier spectrum of 7^. Along side C*{t(T)) is the saturated estimate ttt”^ 
from inequalities and j22| . 
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T = 5 


i = 1.396 


T* -tKi 0.022 


T = 10 


t = 1.41778 
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E{k) 
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t = 1.418002734785 
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Fig. 3. Snapshots of spectra for r = 5, 10 and 25 (top left and right and bottom left respectively), on a lin-log scale. These first three figures show two resolutions 
{N = 1024 and 4096), with and without hyperviscosity, demonstrating the need to control floating point round-off error at small scales, r = 5 shows the full spectrum, 
including dealiased filtered modes to show the small scale error. Thereafter, r = 10 and 25 plots show only the first 500 modes to make clear the initial burst and the onset 
of the slow Gaussian spectrum. The final frame (bottom right) shows only the N = 1024 case with hyperviscosity, now with curves at r = 100, 200 and 500 showing the 
slow broadening of the spectrum at late times. 


terms of the spectrum: 


ll7(-,t)l|co < 


Sk 


= G ZnHri U Zf, 


TT k, 


^odd ' 
k —1 OO. 


[19] 

k-ll2< |k| < fc + 1/2} 


by combining the asymptotic formulae (141 and (15l: 


(7(-,i) ) 


2r 


Il7(-,^)llc 




r.(t) + l 




ll7(',i)l|oo < y/nC{t)S{ty 


[ 20 ] 


The special condition on odd-odd or even-even modes is due 
to the discrete sym metr y of o ur i nitial condition. 

Replacing the fit ( |18[ ) into (191 leads to a bound involving 
an infinite sum over k with an ultraviolet divergence in the 
limit of small 5{t). We can approximate this as follows [12| : 


The left-hand-side of this equation can be written in terms of 
the energy spectrum, s o if we follow similar steps as in the 
derivation of inequality (201 we obtain (7(-,t)^) ~ ^=^<5(1)“^. 
Therefore we get 


C{t): 




[ 21 ] 


Using this we can go back to inequality (201 and show that it 
is equivalent to: 

^ [22] 


where T is the gamma (factorial) function. A further improve¬ 
ment is obtained by noticing the behaviour of n{t) at late times 
from figure 1^ where it is clear that n{t) —>■ —1. Therefore we 
obtain, in tms limit, 




This inequality is in fact saturated, as confirmed by our nu¬ 
merical simulation (figurej^ . In terms of original time variable 
we get 


m< 


This inequality alone cannot be used to estimate the behaviour 
of S{t), since the independent factor C{t) is involved as well, so 
an extra equation is needed. This extra equation is provided 


-W-ii-TT{T* -t)) y -\n{n{T* -t))’ 


which illustrates that the loss of regularity is very slow in 
the original time variable. In fact if one were to consider 
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the reliability time with the saturated spectrum one would 
find that Trei ~ so that for N = 256, Trei « 2 x 10® or 
T* 10“®°°°. 


Errors. Using the definition in the ‘normalised’ L® norm of 
the error (not relative error) is given by 


Q{f,9) 


ll/-g||2 

ll/lb + ||g||2- 


We consider the error associated with the local quan¬ 
tity ||7(-U)l|oo and the global ones (7map), (7^) via = 

Qdl ^nuni ('; II oo 5 ||7ana(', t) || oo ) and Q ^^2 Q{(7num) 5 (7ana) )! 

etc. where the subscripts “num” and “ana” stand for “numer¬ 
ical” and “analytic”. 

The numerical solution of the mapped system does not pro¬ 
vide direct access to the original variable |j 7 (- , t)||oo so the 
following expression is required [5] 


Il7(-,t(r))|| 


Il7o||ooexp 



[23] 


where /Q^( 7 map)dr' is computed using Simpson’s rule. We 
compare both this mapped estimate and the direct supremum 
norm from the original system against the analytical solution 
of Eq. ( |11[ ). Care is taken in solving Eq. so that the 
time steps from the original system are used to solve on in¬ 
tervals which coincide with the data points and arbitrary pre¬ 
cision of the required level is used. As shown in the previ¬ 
ous section the solution remains well resolved spatially, even 
at relatively modest resolutions, therefore we omit the error 
study of spatial convergence here. We do however show con¬ 
vergence with respect to timestep dr in figure at resolution 
N = 1024. Overall we observe an exponential growth (in r) 



Fig. 4. Time evolution of the error measures Q-y, (top) Q(-y2 ^ and Q(^^ 2 ^ 

(bottom) showing the convergence with time step dr. 


of Q-y from the original system, compared to an almost uni- 
form error (~ 10“^°) from the mapped version at converged 
dr. Interestingly convergence occurs in the mapped system at 
a smaller level of dr. 

In contrast to the earlier result in [^, we find that both Qj 
and ^ behave similarly; recall that (7map) is the primary 

variable for the assessment of the local quantity ||7(-,t)||oo in 
the mapped system (equation (23l). This implies that in this 
case (where (7®) is not an invariant), Qj is simply slaved to 
Q(-y 2 y Ref. [ 5 ] contains a detailed discussion on error sources 
in the mapped and original variables, and highlights some sub¬ 
tleties surrounding the behaviour of {7®). Here the situation 
is somewhat more straightforward: the original system con¬ 
tains unbounded error growth due to a fundamental loss of 
precision in the independent variable. This is explained by 
the earlier convergence and the ‘saturation’ of error near the 
double-precision limit (r ~ 37). 


Assessing blowup time T* and proximity to it, T* — t. Pre¬ 
vious methods for assessing the value of T*, e.g. fitting the 
behaviour of || 7 (' , t)||oo to a power law (T* —1)“, are based on 
the assumption that the solution is incurring significant errors 
at intermediate time and T* requires careful extrapolation. 
Here the solution is remaining well resolved until late times 
and we find that using the original system with the adaptive 
timestep given above, t converges to T* to within ~ 10“®’^. 
This accuracy is surprising given it arises from a simple sum 
ti = '^''"/Il7('i ^'^d it cannot be improved by fitting 

or even by arbitrary precision arithmetic to sum dt which are 
below the machine precision threshold. Understanding this 
accuracy is aided by attempting the comparable exercise for 
the mapped system. Here the recovery of t is given by 


t(r) = 


f 


II7('U(t 


-dr'. 


where ||7(-,t(r'))||oo is obtained from formula (23l. This inte¬ 
gral should converge to T* as r — >■ 00. Numerically computing 
it results in a saturation of error ~ 10“° (slaved to the error 
in ||7(-,r)||oo) when dr = 10“^^ for sufficiently large r. This 
almost leads to a paradox: why should a quantity with lower 
late time error produce a poorer estimate for T* when the 
procedure for the estimate is qualitatively the same. The rea¬ 
son is that it is the early errors which pollute the estimate for 


142 



Fig . 5. Relative error of proximity P(t) compared to the asymptotic formula 
jl3| plotted against T* —t (bottom axis) and r (top axis). Curves show convergence 
in dr and an error rvj 10~^ persistent to exceptionally small values of T* — t. 
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t(r) as these are larger in the mapped system and occur at a 
point where they will contribute more significantly to the final 
integral (dr/|| 7 (- , l)||oo is largest). 

However, one should proceed with caution when dismissing 
the ability of the mapped system at assessing its original tem¬ 
poral position: error in the assessment of T* is not to be con¬ 
fused with error in the proximity to T*. Although the original 
system can integrate to within lO'^'^ of T*, it is impossible 
to assess any behaviour beyond this point: this is a solid bar¬ 
rier for the method. For the mapped system this is not the 
case: through t(r) we can produce an estimate for T* —t as a, 
function of t by considering the following ‘proximity’ integral 


T* -t 



1 


dr'. 


Recalling that ||7(-,t)||oo is recovered from (7map) in the 
mapped system via formula ( |23| ), we fit the behaviour of this 
global measure in the preceding (r —10) window via the ansatz 
(7map) ~ ~ Inserting this into what is now the double 

integral for T* — t we obtain 


(1-2k)t 2m /I _ o \2m 

Il 7 ('l ^)l|o<D 


slow spatial collapse and unusual (Gaussian) Fourier spectrum 
very near singularity time. This means that the solution will 
remain well spectrally converged until extraordinarily close to 
singularity time for even modest resolutions. In turn this im¬ 
plies a fundamental constraint on the original system: in the 
original variables one can only hope to approach T* to the 
precision of the floating point arithmetic being used, usually 
double-precision, ~ 10“^®. Because of this lack of digits in the 
independent temporal variable, assessing any quantities from 
the original system is a hazardous undertaking as errors grow 
exponentially. In other words, not only does the proximity to 
T* present a floating point barrier: it also harms the accurate 
assessment of the late time behaviour of the system before the 
barrier is reached. On the other hand, the mapped system has 
no floating point arithmetic barrier as the singularity time is 
now at infinity and we observe uniform errors until T* — t is 
exceptionally small in the figures shown). 


Another floating point barrier also becomes apparent, 
namely that ||7||oo will eventually overflow, i.e. exceed ~ 10^°® 
at r « 715. Luckily the mapping allows us to postpone this 
barrier further by simply computing log ||7||oo ( i.e. , outputting 
the exponent of the right hand side of equation (23 l)and use an 
arbitrary precision exponential in post-processing if required. 


where F is the incomplete gamma function. This provides a 
running estimate for T* —t w hich we can validate against the 
asymptotic formula, equation (131. Figurel^shows the relative 
error in P{t) as a function of f in oiRer to demonstrate 
how the error depends on the absolute proximity to T*. We 
find relative errors of the order of 10~^ persisting far beyond 
the machine precision limit, and converging at larger dr than 
Q-y, presumably due to the accuracy of the fitting procedure. 


There is already some evidence that, depending on the type 
of initial conditions, the full 3D problem has a changing late¬ 
time regime where either a depletion of nonlinearity slows vor- 
ticity growth m or the collision of two vortex sheets accel¬ 
erates the loss of regularity m- It is therefore hoped that 
mapping the full 3D problem will give renewed confidence in 
the late-time behaviour of the next generation of simulations 
of 3D Euler. 


Conclusion and Discussion 

In this p^er we have shown that only by mapping the singular 
system to a regular one ([^, Q, can certain uncon¬ 

ventional late-time behaviours be observed and asymptotic 
trends be established. The first unusual feature shown is the 
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